Astronomical pacing of the Lau biogeochemical Event captured in the Kosov Quarry section
A study on the imprint of astronomical cycles in the Kosov quarry spanning the Lau Event during the Silurian
Authors
Michiel Arts1
Anne-Christine da Silva1
1 - SEDICLIM lab, Department of Geology, University of Liege, Belgium
Published
October 1, 2026
Note
This study is based onthe induration record extracted from the litholog of Fryda, Jiri & Lehnert, Oliver & Frýdová, Barbora & Farkas, Juraj & Kubajko, Michal. (2020). Carbon and sulfur cycling during the mid-Ludfordian anomaly and the linkage with the late Silurian Lau/Kozlowskii Bioevent. Palaeogeography, Palaeoclimatology, Palaeoecology. 564. doi:10.1016/j.palaeo.2020.110152
The Kosov Quarry section in the Prague Basin preserves one of the most complete Silurian successions spanning the Lau biogeochemical Event. This interval records major perturbations in the carbon cycle, climate cooling, redox fluctuations, and faunal turnover, yet the pacing and duration of these changes remain poorly constrained. Leveraging its exceptional completeness, we developed a new astrochronology based on the induration record, which allows the identification of cycles ranging from the half-precession to the 405-kyr eccentricity cycle. The astrochronology indicates that the Lau biogeochemical Event started at 424.06 ± 0.55 Ma and lasted 0.87 ± 0.03 Myr, while the associated Lau-Kozlowski Bio-event was much shorter at 30 ± 10 kyr. Astronomical forcing is expressed not only in the induration record but also in the rate of change (‰/kyr) of the δ13Ccarb curve and in the lag-1 sea-level proxy. The δ13Ccarb rate-of-change record shows a strong 100-kyr eccentricity imprint during the onset of the Lau Event, reaching a maximum of 0.09 ‰/kyr. Similarly, the lag-1 sea-level record is dominated by the 405-kyr eccentricity cycle and appears in-phase with insolation, consistent with (405-kyr) eccentricity-paced glacio-eustatic sea-level fluctuations during the Ludlow.
1. Introduction
The Lau biogeochemical Event is associated with one of the largest carbon cycle perturbations of the entire Phanerozoic. The Lau biogeochemical Event is marked by a large positive δ¹³C excursion, faunal turnover, and evidence for sea-level and redox fluctuations (Frýda and Manda, 2013; Frýda et al., 2020). These changes indicate a tight coupling between the carbon cycle, climate, and oceanographic processes; however, the tempo, pacing, and internal structure of the Event remain poorly constrained. Without a precise temporal framework, it is not possible to determine whether the Lau reflects a brief perturbation or a protracted interval of instability.
To address this, we will conduct a cyclostratigraphic study on the Kosov quarry section, which provides a nearly complete archive, allowing us to establish a high-resolution astrochronological framework for the Lau biogeochemical Event. The induration record will form the basis for constructing the astrochronology, which will enable us to assign durations and associated uncertainties to the Lau biogeochemical Event and its associated (sub)intervals. The lag-1 sea-level proxy will be derived from the tuned induration record, which will provide constraints and insights into sea-level variability during the Lau biogeochemical event. The rate of change (‰/kyr) calculated from the δ13Ccarb record will be used to gain insights into carbon cycle dynamics during the Lau biogeochemical Event. Furthermore, the imprint of astronomical cycles on the δ13Ccarb rate of change (‰/kyr) and lag-1 sea-level will be delineated to understand the role that astronomical cycles might have played in pacing climate, organic carbon burial, and sea-level variations during the Lau biogeochemical Event.
1.1 Geological setting
Our study focuses on upper Silurian (Ludfordian) strata of the Prague Basin, a sedimentary succession within the Barrandian (Teplá-Barrandian terrane) (Frýda and Manda, 2013; Frýda et al., 2020, 2021) (Figure 1). During the Ludfordian, this block was located at ~25–30°S (Tasáryová et al., 2014; Scotese, 2021). Facies distribution indicates hemipelagic deposition surrounded by shallow- to deeper-marine environments (HORNÝ, 1955; Kříž, 1991). It remains debated whether the Barrandian’s geodynamic setting belongs to the peri-Gondwanan terrains or is its own isolated Perunica microcontinent (Fatka and Mergl, 2009).
Code
par(mar =c(4, 2, 2, 2))# Empty plot to define margins and space in the panelplot(1, type ="n", axes =FALSE, xlab ="", ylab ="", main="A")# Get the coordinates of the current plot regionusr <-par("usr") # c(x1, x2, y1, y2)# Draw the image to fill the entire panelrasterImage(img_2, usr[1], usr[3], usr[2], usr[4])
Figure 1. Regional overview of the Kosov section. Copied from Fryda et al. (2020)
Code
par(mar =c(4, 6, 2, 2))
The studied Kosov section (No.JF195; 49°56′12.2″N, 14°03′20.1″E) provides one of the most complete global records of the mid-Ludfordian carbon isotope excursion (Frýda and Manda, 2013; Frýda et al., 2020, 2021) (Figure 2).The succession comprises carbonates with marly interbeds and shales, including micritic mudstones, skeletal limestones, and crinoidal grainstones (Kříž, 1992; Lehnert et al., 2007). The site has been widely studied for its sedimentology, palaeontology, and faunal changes across the Lau/Kozlowskii Bioevent (LKB) and Mid-Ludfordian Carbon Isotope Excursion (MLCIE) (Kříž, 1992; Štorch, 1995; Manda and Kříž, 2006; Frýda and Manda, 2013).
Figure 2. The original litholog, chronostratigraphy, event zonation, induration record and the δ13Ccarb curve. A. The Litholog of Frýda et al. (2020). B. Chronostratigraphic subdivision. C. Event stratigraphic intervals. D. The induration record extracted from the litholog E. The δ13Ccarb record of Frýda et al. (2020)
2. Materials and Methods
A cyclostratigraphic analysis will be conducted on the induration record (see SI for the R code). The analyses will leverage the function of the WaverideR (Arts, 2023; Arts et al., 2024) astrochron R packages (Meyers, 2015, 2019) to conduct the cyclostratigraphic study.
The cyclostratigraphic analysis begins by applying the Time Optimisation (TimeOpt) function of the astrochron R package to the induration record, yielding a first-order estimate of the sedimentation rate (Meyers, 2015). The TimeOpt function will be run twice, once to evaluate the 405-kyr eccentricity amplitude modulation of the 100-kyr eccentricity band, and the concentration of power at the 100-kyr and 405-kyr eccentricity frequencies, and once to evaluate the eccentricity amplitude modulation of the precession band, and the concentration of power at the precession and eccentricity frequencies. The optimal fit of the astronomical cycles will be evaluated for a range of sediment accumulation rates spanning between 0.1 and 4 cm/kyr.
The first-order sedimentation rate estimate from the TimeOpt runs is used to extract the 405-kyr eccentricity cycle from the induration record. The extracted cycle is then used to conduct a minimal tuning, where the distance between two successive peaks of the extracted cycle is set to the duration of the interpreted astronomical cycle. This duration is subsequently converted into an age model and a sedimentation-rate curve. The sedimentation-rate curve is then used to convert cm/kyr to the period (m) of different astronomical cycles. These cycles will then be plotted on top of the CWT scalogram to verify whether the curves pass through areas of high spectral power, validating the presence of astronomical cycles. If this is the case, the period (m) of the 405-kyr eccentricity cycle will be tracked in the CWT.
The tracked period (m) curve of the 405-kyr eccentricity cycle can then be used to create a floating astrochronology with an assigned uncertainty based on the analytical uncertainty of the Wavelet (Arts et al., 2024). This floating age model with uncertainty is then anchored to the astrochronologically calibrated age of the Ludlow-Pridoli boundary (423.03 ±0.53 Ma) (Arts et al., 2024). The anchored numerical age model can then be used to assign ages and durations (with uncertainty) to the record of the Kosov quarry section and its subsections.
The mean age model is then used to convert the induration and δ13Ccarb record to the time domain. Allowing for the study of the imprint of astronomical cycles in the induration record in the time domain. Special attention will be paid to the phase relationship between the 100-kyr eccentricity cycle directly extracted from the record and the 100-kyr eccentricity cycle extracted from the Hilbert transform of the precession cycle, and the phase relationship between the 405-kyr eccentricity cycle directly extracted from the record and the 405-kyr eccentricity cycle extracted from the Hilbert transform of the 100-kyr eccentricity cycle. The phase relationship between these cycles will allow us to infer the phase relationship between cycles extracted from the proxy record and those the true phase of the astronomical cycle, as well as its relationship with changes in insolation (Hinnov, 2000; Laskar et al., 2004).
The rate of change (‰/kyr) will be calculated from the δ13Ccarb record, which has been converted to the time domain. Which again will be investigated for the imprint of astronomical cycles. A lag-1 curve was generated based on the lag-1 autocorrelation coefficient using a windowed analysis Monte-Carlo on the induration record in the time domain. This record will function as a proxy record for the sea-level curve (Li et al., 2018). The imprint of astronomical cycles in the lag-1 record will be investigated and compared to the imprint of astronomical cycles in the induration recordd.
3. Results
3.1. TimeOpt-derived sedimentation rates
The TimeOpt function was run twice, once to optimise the 405-kyr eccentricity amplitude modulation of the 100-kyr eccentricity band, and the concentration of power at the 100-kyr and 405-kyr eccentricity frequencies and once to optimise the eccentricity amplitude modulation of the precession band, and the concentration of power at the precession and eccentricity frequencies. The TimeOpt results indicate that the optimal sedimentation rate is around 3 cm/ka, with the peak spanning sedimentation rates between 2 and 4 cm/ka (see Figure 3).
Code
#run timeOpt#we need to run it twice for a correct plot later ontargetP_AM=c(405.7, 130.7, 123.8, 98.9,94.9)targetP=c(21,17)res1=timeOpt(dat = induration,sedmin =0.1,sedmax =4,numsed =500,fit =2,fitModPwr =TRUE,targetE=targetP_AM,roll =10^4,output=1)res2 <-timeOpt(dat = induration,sedmin =0.1,sedmax =4,numsed =500,fit =2,fitModPwr =TRUE,targetE=targetP_AM,roll =10^4,output=2,genplot=TRUE)res3=timeOpt(dat = induration,sedmin =0.1,sedmax =4,numsed =500,fit =1,fitModPwr =TRUE,targetP=targetP,roll =10^4,output=1)res4 <-timeOpt(dat = induration,sedmin =0.1,sedmax =4,numsed =500,fit =1,fitModPwr =TRUE,targetP=targetP,roll =10^4,output=2,genplot=TRUE)
Code
layout(matrix(c(1, 2, 3, 4,5,6,7,8), nrow =4, ncol =2, byrow =TRUE)) par(mar =c(4, 4, 3, 3))# --- Panel A ---plot(res1[, 1], res1[, 2], cex =0.75, cex.lab =1.2, cex.main =1.3, col ="red", xlab ="", ylab ="", main =expression(paste(bold("Fit: "), {"r"^2 }["envelope"], " (red) and ", {"r"^2 }["power"], " (gray)")))axis(2, col.axis ="red")mtext(expression("r"^2), side =2, line =2, cex =0.9)mtext("Sedimentation rate (cm/ka)", side =1, line =2.3, cex =0.8)par(new =TRUE)plot(res1[, 1], res1[, 3], col ="#00000064", xlab ="", ylab ="", type ="l", axes =FALSE, lwd =2)axis(4, ylim =c(0, max(res1[, 3])), lwd =1, col ="black")usr <-par("usr")text(usr[1], usr[4] -0.05*(usr[4]-usr[3]), "A", pos =4, font =2, cex =1.5)# --- Panel B ---plot(res1[, 1], res1[, 4], type ="l", lwd =2, cex.lab =1.2, cex.main =1.3, col ="black", xlab ="", ylab ="", main =expression(paste(bold("Optimal Fit: "), {"r"^2 }["opt"])))mtext(expression({"r"^2}["opt"]), side =2, line =1.9, cex =0.9)mtext("Sedimentation rate (cm/ka)", side =1, line =2.3, cex =0.8)usr <-par("usr")text(usr[1], usr[4] -0.05*(usr[4]-usr[3]), "B", pos =4, font =2, cex =1.5)# --- Panel C ---ylim1 =c(min(res2[,4], res2[,5]), max(res2[,4], res2[,5]))plot(res2[,1], res2[,4], cex =0.5, cex.lab =1.2, cex.main =1.2, xlab ="", ylab ="", main ="Envelope (red); Reconstructed Ecc. Model (black)", col ="red", type ="l", ylim = ylim1)lines(res2[, 1], res2[,5], lwd =1.5)mtext("Std. Value", side =2, line =2, cex =0.9)mtext("Time (ka)", side =1, line =2.3, cex =0.8)usr <-par("usr")text(usr[1], usr[4] -0.05*(usr[4]-usr[3]), "C", pos =4, font =2, cex =1.5)# --- Panel D ---ylim1 =c(min(res2[, 3], -1* res2[, 4]), max(res2[, 3], res2[, 4]))plot(res2[,1],res2[,3], col ="blue", cex =0.5, cex.lab =1.2, cex.main =1.3, xlab ="", ylab ="", main ="Envelope (red); Filtered Data (blue)", ylim = ylim1)lines(res2[,1],res2[,3], col ="blue")lines(res2[,1],res2[,4], col ="red")lines(res2[,1], -1*res2[,4], col ="red")abline(h =0, col ="black", lty =3)mtext("Std. Value", side =2, line =2, cex =0.9)mtext("Time (ka)", side =1, line =2.3, cex =0.8)usr <-par("usr")text(usr[1], usr[4] -0.05*(usr[4]-usr[3]), "D", pos =4, font =2, cex =1.5)# --- Panel E ---plot( res3[, 1], res3[, 2], cex =0.75, cex.lab =1.2, cex.main =1.3, col ="red", xlab ="", ylab ="", main =expression(paste(bold("Fit: "), {"r"^2 }["envelope"], " (red) and ", {"r"^2 }["power"], " (gray)")))axis(2, col.axis ="red")mtext(expression("r"^2), side =2, line =2, cex =0.9)mtext("Sedimentation rate (cm/ka)", side =1, line =2.3, cex =0.8)par(new =TRUE)plot( res3[, 1], res3[, 3], col ="#00000064", xlab ="", ylab ="", type ="l", axes =FALSE, lwd =2)axis(4, ylim =c(0, max( res3[, 3])), lwd =1, col ="black")usr <-par("usr")text(usr[1], usr[4] -0.05*(usr[4]-usr[3]), "E", pos =4, font =2, cex =1.5)# --- Panel F ---plot( res3[, 1], res3[, 4], type ="l", lwd =2, cex.lab =1.2, cex.main =1.3, col ="black", xlab ="", ylab ="", main =expression(paste(bold("Optimal Fit: "), {"r"^2 }["opt"])))mtext(expression({"r"^2}["opt"]), side =2, line =1.9, cex =0.9)mtext("Sedimentation rate (cm/ka)", side =1, line =2.3, cex =0.8)usr <-par("usr")text(usr[1], usr[4] -0.05*(usr[4]-usr[3]), "F", pos =4, font =2, cex =1.5)# --- Panel G ---ylim1 =c(min( res4[,4], res4[,5]), max( res4[,4], res4[,5]))plot( res4[,1], res4[,4], cex =0.5, cex.lab =1.2, cex.main =1.2, xlab ="", ylab ="", main ="Envelope (red); Reconstructed Ecc. Model (black)", col ="red", type ="l", ylim = ylim1)lines( res4[, 1], res4[,5], lwd =1.5)mtext("Std. Value", side =2, line =2, cex =0.9)mtext("Time (ka)", side =1, line =2.3, cex =0.8)usr <-par("usr")text(usr[1], usr[4] -0.05*(usr[4]-usr[3]), "G", pos =4, font =2, cex =1.5)# --- Panel H ---ylim1 =c(min( res4[, 3], -1* res4[, 4]), max( res4[, 3], res4[, 4]))plot( res4[,1], res4[,3], col ="blue", cex =0.5, cex.lab =1.2, cex.main =1.3, xlab ="", ylab ="", main ="Envelope (red); Filtered Data (blue)", ylim = ylim1)lines( res4[,1], res4[,3], col ="blue")lines( res4[,1], res4[,4], col ="red")lines( res4[,1], -1* res4[,4], col ="red")abline(h =0, col ="black", lty =3)mtext("Std. Value", side =2, line =2, cex =0.9)mtext("Time (ka)", side =1, line =2.3, cex =0.8)usr <-par("usr")text(usr[1], usr[4] -0.05*(usr[4]-usr[3]), "H", pos =4, font =2, cex =1.5)
Figure 3. TimeOpt results. A. Fit: r²[envelope] (red) and r²[power] (grey) for TimeOpt run 405 vs 100-kyr ecc. B. Optimal Fit: r²[opt] for TimeOpt run 405 vs 100-kyr ecc. C. Envelope (red); Reconstructed Ecc. Model (black) for TimeOpt run 405 vs 100-kyr ecc. D. Envelope (red); Filtered Data (blue) for TimeOpt run 405 vs 100-kyr ecc. E. Fit: r²[envelope] (red) and r²[power] (grey) for TimeOpt run 100-kyr eccentricity versus precession. F. Optimal Fit: r²[opt] for TimeOpt run 100-kyr eccentricity versus precession. G. Envelope (red); Reconstructed Ecc. Model (black) for TimeOpt run 100-kyr ecc. versus precession. H. Envelope (red); Filtered Data (blue) for TimeOpt run 100-kyr eccentricity versus precession.
3.2 Extract the 405-kyr cycle and generate a sedimentation-rate curve
Based on the TimeOpt results, the 405-kyr eccentricity cycle is extracted from the induration record using a Taner bandpass filter, which spans from 4.86 to 16.2, corresponding to the peak of the optimal sedimentation rate range (1.2-4 cm/ka) as determined by the TimeOpt analysis results. Next, a sedimentation-rate curve is constructed using the minimal tuning technique, in which the distance between two successive peaks of the extracted cycle is set to the duration of the interpreted astronomical cycle. This duration is subsequently used to generate the sedimentation-rate curve (see Figure 4).
Figure 4. Minimal tuning. A. Induration record and the extracted 405-kyr eccentricity cycle. B. Sedimentation rate and sedimentation rate (cm/kyr) based on the minimal tuning.
3.3 Tracking the period (m) of the 405-kyr eccentricity cycle in the the CWT scalogram
The sedimentation rate curve attained from the minimal tuning was then used to convert cm/kyr to period (m) of different astronomical cycles. These cycles were then plotted on top of the CWT scalogram (see figure 4). We can see that the curves plot through areas of high spectral power corresponding to astronomical cycles that should be present according to the TimeOpt-derived sedimentation rate modelling runs. Next, the period (m) of the 405-kyr eccentricity cycle was tracked in the CWT.
Figure 5. The Continuous Wavelet Transform (CWT) of the induration record extracted from the litholog of Frýda et al. (2020) and the period (m) of astronomical cycles based on the sedimentation rates attained from the minimal tuning and the tracked period (m) of the 405-kyr eccentricity cycle. Black dotted = tracked period (m) the 405-kyr eccentricity cycle Red = the 405-kyr eccentricity cycle (minimal tuning based) Purple = the 100-kyr eccentricity cycle (minimal tuning based) Blue = the obliquity cycle (minimal tuning based) Dark green = the precession cycle (minimal tuning based)
3.4 Combine the tracked period (m) and a tie-point curve to create an anchored astrochronoloy
The tracked period (m) curve of the 405-kyr eccentricity cycle was then be used to create a floating astrochronology with an assigned uncertainty based on the analytical uncertainty of the Wavelet (Arts et al., 2024). This floating age model with uncertainty was then anchored to the astrochronologically calibrated age of the Ludlow-Pridoli boundary (423.03 ±0.53 Ma) (Arts et al., 2024). This anchored numerical age model could then be used to assign ages and durations (with uncertainty) to the record of the Kosov quarry section and its subsections
Figure 6. Astrochronological numerival age- depth model.The black line is the mean age-depth curve, and the red and blue lines are time plus and minus two standard deviations (2σ)
3.5 Ages and durations for intervals including uncertainty
The age model provides durations with uncertainty for the chrono, litho-, and event-stratigraphic intervals previously identified in the succession (Frýda and Manda, 2013; Frýda et al., 2021). The calculated ages, durations, and their associated uncertainties were rounded to the nearest 10 kyr to align with the precision of the anchor point.
Code
kosov_time<-cbind(curve2time_unc_anchor_res[[1]],rowMeans(results))kosov_sd<-rowSds(as.matrix(results))for(i in1:nrow(my_data)){row_nr_1 <- DescTools::Closest(kosov_time[,1],my_data[i,2],which =TRUE)row_nr_2 <- DescTools::Closest(kosov_time[,1],my_data[i,3],which =TRUE)my_data[i,4] <-round(mean(as.matrix(results[row_nr_1,]),na.rm =TRUE)/10,0)/100my_data[i,5] <-round(2*sd(as.matrix(results[row_nr_1,]),na.rm =TRUE)/10,0)*10my_data[i,6]<-round(mean(as.matrix(results[row_nr_2,]),na.rm =TRUE)/10,0)/100my_data[i,7] <-round(2*sd(as.matrix(results[row_nr_2,]),na.rm =TRUE)/10,0)*10my_data[i,8] <-round(mean(as.matrix(na.omit(results[row_nr_1,]-results[row_nr_2,])))/10,0)*10my_data[i,9] <-round(2*sd(results[row_nr_2,]-results[row_nr_1,],na.rm =TRUE)/10,0)*10}for(i in1:nrow(my_data)){for( j in4:ncol(my_data)){if (my_data[i,j]<=10){ my_data[i,j] <-10 }}}bottom_record <- kosov_time[DescTools::Closest(kosov_time[,1],my_data[1,2],which=TRUE),2]base_pridoli <- kosov_time[DescTools::Closest(kosov_time[,1],my_data[8,2],which=TRUE),2]Top_record <- kosov_time[DescTools::Closest(kosov_time[,1],my_data[1,3],which=TRUE),2]base_LAU<- kosov_time[DescTools::Closest(kosov_time[,1],my_data[4,2],which=TRUE),2]base_S_zone<- kosov_time[DescTools::Closest(kosov_time[,1],my_data[5,2],which=TRUE),2]base_F_zone<- kosov_time[DescTools::Closest(kosov_time[,1],my_data[6,2],which=TRUE),2]Top_F_zone<- kosov_time[DescTools::Closest(kosov_time[,1],my_data[6,3],which=TRUE),2]base_event <- kosov_time[DescTools::Closest(kosov_time[,1],my_data[2,2],which=TRUE),2]# Create gt table with header and custom column label (no pipe)tab <-gt(my_data)tab <-cols_label(tab,position_bottom ="position (m) bottom",position_top ="position (m) top",age_bottom ="Age bottom (Ma)",uncertainty_age_bottom="uncertainty bottom (kyr) +/-2sd ",age_top ="Age top (Ma)",uncertainty_age_top="uncertainty top (kyr) +/-2sd ",duration ="duration (kyr)",uncertainty ="uncertainty in duration (kyr) +/-2sd ")tab <-tab_header(tab, title ="Table 1. Ages and durations based on the anchored astrochronological age model")tab
Table 1. Ages and durations based on the anchored astrochronological age model
Interval
position (m) bottom
position (m) top
Age bottom (Ma)
uncertainty bottom (kyr) +/-2sd
Age top (Ma)
uncertainty top (kyr) +/-2sd
duration (kyr)
uncertainty in duration (kyr) +/-2sd
Entire record
-7.2
25.0
424.60
560
422.80
530
1800
190
Lau biogeochemical Event
-0.4
20.0
424.06
550
423.19
530
870
90
LKB Event (trace metal peak)
-0.4
0.0
424.06
550
424.03
550
30
10
R-zone
0.0
1.4
424.03
550
423.94
540
90
10
S-Zone
1.4
11.6
423.94
540
423.64
540
310
30
F-Zone
11.6
20.0
423.64
540
423.19
530
440
50
Ludlow part section
0.0
22.0
424.03
550
423.03
530
1000
110
Pridoli part section
22.0
25.0
423.03
530
422.80
530
230
20
3.6 Extract and plot the record and the astronomical cycles in the time domain
The age model was used to convert the induration and δ13Ccarb record to the (numerical) time domain. Next, the CWT was performed on the tuning induration record. Peaks corresponding to the 8.5-kyr half precession, ~16-kyr precession, 28-kyr bliquity,50-kyr eccentricity, 100-kyr eccentricity, 200-kyr eccentricity, 405-kyr eccentricity and the 1.6-Myr eccentricity cycle can be observed in the CWT scalogram. The precession, obliquity, 100-kyr eccentricity and 405-kyr eccentricity were also extracted from the record. The Hilbert transform was applied to the precession, obliquity and 100-kyr eccentricity cycles. The 405-kyr eccentricity cycle was also extracted from the Hilbert transform-derived amplitude record of the 100-kyr eccentricity and compared to the 405-kyr eccentricity directly extracted from the record, showing an anti-phased relationship.
Figure 7. Induration proxy record in the time domain, including the wavelet scalogram of the induration record and astronomical cycles extracted from said record. A Stages. B. Event zone subdivision C. The δ13Ccarb record. D. The Induration record. E. Wavelet scalogram of the induration record with the average spectral power on top. The black vertical lines in the wavelet scalograms are durations of known astronomical cycles. From left to right, these cycles are the 8.5-kyr half precession, ~16-kyr precession, 28-kyr obliquity,50-kyr eccentricity, 100-kyr eccentricity, 200-kyr eccentricity, 405-kyr eccentricity and the 1.6-Myr eccentricity cycle. F. The black line is a 405-kyr eccentricity cycle extracted from the Induration record. The red line is the 405-kyr eccentricity cycle extracted from the Hilbert transform of the 100-kyr eccentricity cycle of the Induration record. G. The 100-kyr eccentricity cycle was extracted from the Induration record (black line) and the Hilbert transform of the 100-kyr eccentricity cycle (red line). H.) Obliquity cycle extracted from the Induration record (black line), and the Hilbert transform of the obliquity cycle (red line). I. The precession cycle extracted from the Induration record (black line) and the Hilbert transform of the precession cycle (red line).
3.7 The phase relationships
The 405-kyr eccentricity cycle was extracted from the induration record and extracted from the Hilbert transform of the 100-kyr eccentricity cycle of the induration record. The 100-kyr eccentricity cycle was extracted from the induration record and extracted from the Hilbert transform of the precession cycle of the induration record. The relationship between the 405-kyr eccentricity cycle directly extracted from the record and the one obtained from the Hilbert transform of the 100-kyr cycle shows an anti-phased relationship. The relationship between the 100-kyr eccentricity cycle directly extracted from the induration record and that obtained from the Hilbert transform of the precession cycle demonstrates a generally anti-phased relationship as well.
Code
ylims <-c(min(kosov_time[, 2])-5,max(kosov_time[,2])+5)ylims <-rev(ylims/1000)layout.matrix <-matrix(c(1,2),nrow =1,ncol =2 ,byrow =TRUE)graphics::layout(mat = layout.matrix,heights =c(1,1),# Heights of the two rowswidths =c(1,1))par(mar =c(4, 4, 4, 2))plot(x = induration_filt_time_405[, 2],y = induration_filt_time_405[, 1]/1000,type ="l",ylim =ylims,xlab ="405-kyr ecc",col="black",ylab ="Age (Ma)",main="A")lines(x = (induration_filt_time_110_hilb_405[, 2]-mean(induration_filt_time_110_hilb_405[, 2]))*2+mean(induration_filt_time_405[,2]), y = induration_filt_time_110_hilb_405[, 1]/1000, col ="red")plot(x = induration_filt_time_110[, 2],y = induration_filt_time_110[, 1]/1000,type ="l",ylim =ylims,xlab ="100-kyr ecc",col="black",ylab ="Age (Ma)",main="B")lines(x = (induration_filt_time_prec_hilb_110[, 2]-mean(induration_filt_time_prec_hilb_110[, 2]))*2+mean(induration_filt_time_110[,2]), y = induration_filt_time_prec_hilb_110[, 1]/1000, col ="red")
Figure 8. Astronomical cycles extracted from the tuned induration record. A. The black line is a 405-kyr eccentricity cycle extracted from the Induration record. The red line is the 405-kyr eccentricity cycle extracted from the Hilbert transform of the 100-kyr eccentricity cycle of the Induration record. B. The black line is a 100-kyr eccentricity cycle extracted from the Induration record. The red line is the 100-kyr eccentricity cycle extracted from the Hilbert transform of the precession cycle of the Induration record.
3.8 Rate of change in the δ13Ccarb curve
The rate of change (‰/kyr) of the δ13Ccarb curve was calculated from the LOWESS regression applied to the δ¹³Ccarb record. The regression was applied to minimise the influence of short-term local noise (Figure 9). The rate of change (‰/kyr) reaches a maximum of 0.09 (‰/kyr) during the onset of the Lau biogeochemical Event. The rate of change (‰/kyr) record contains a strong imprint of the 100-kyr eccentricity cycle, which was subsequently extracted from the record and plotted next to the 100-kyr eccentricity cycle extracted from the induration record. During the Lau, the rate of change (‰/kyr) recorded its largest amplitude variations. The 100-kyr eccentricity cycle extracted from the rate of change (‰/kyr) record is generally anti-phased with the 100-kyr eccentricity cycle extracted from the induration record.
Figure 9. Rate of change (‰/kyr) of the δ13Ccarb record. A Stages. B. Event zone subdivision. C. The Induration record. D. The δ13Ccarb record (light green), and the LOWESS smoothed curve (dark green) E. The δ13Ccarb rate of change record (‰/kyr) and the 100-kyr eccentricity cycle extracted from said record F. The 100-kyr eccentricity cycle extracted from the δ13Ccarb rate of change record (‰/kyr) (dark green) and the 100-kyr eccentricity cycle extracted from the induration record.
3.9 The lag-1 sea-level curve
The lag-1 function was used to calculate the lag-1 autocorrelation coefficient using a windowed analysis Monte-Carlo analysis. The resulting curve serves as a proxy for sea-level changes (Li et al., 2018). The lag-1 does show a few trends (see Figure 10). Before the Lau excursion, the curve is characterised by relatively high values superimposed by some low-amplitude fluctuations. At the onset of the Lau event, a minor curve rise is observed, followed by a large drop towards a minimum during the middle of the Lau δ13Ccarb excursion. Toward the top of the studied succession, the curve returns to higher and more regular values. The lag-1 record contains the imprint of the 405-kyr eccentricity cycle, which was extracted and compared to the 405-kyr eccentricity cycle extracted from the induration record. The 405-kyr eccentricity cycle extracted from the rate lag-1 record is generally anti-phased with the 405-kyr eccentricity cycle extracted from the induration record.
This code takes a long time to run so one can also download the results from github induration_time_lag_1 <- lag_1( data = induration_time_2, n_sim = 1000, run_multicore = TRUE, win_max = 300, win_min = 50, verbose = TRUE ) write.csv(induration_time_lag_1,“induration_time_lag_1.csv”)
Figure 10. Rate of change of the δ13Ccarb record. A Stages. B. Event zone subdivision. C. The Induration record. D. The δ13Ccarb record (light green) and the LOWESS smoothed curve (dark green) E. The lag-1 sea-level curve (red) and the 405-kyr eccentricity cycle extracted from said record (brown). F. The 405-kyr eccentricity cycle extracted from the lag-1 curve (brown) and the 405-kyr eccentricity cycle extracted from the Induration record (blue).
4. Discussion
4.1 The astrochronology and the imprint of astronomical cycles
The anchored astrochronological age model was used to assign ages, durations, and associated uncertainties to chrono- and geochronologic units, lithological units, and events in the Kosov quarry section, as tabulated in Table 1. The Lau biogeochemical Event started at 424.06 ± 0.55 Ma and lasted 0.87 ± 0.03 Myr; in contrast, the LKB Event was as short as 30 ± 10 kyr. The durations and associated uncertainties for the Lau Event and its subdivisions are the first astrochronological constraints available for a near-complete record spanning the Lau Event. The durations and associated uncertainties are therefore the best constrained reference values to date. A previous estimate for the onset of the Lau Event is based on the dating of a hiatus in the Cellon section, spanning 423.83 ± 0.55 Ma to 424.19 ± 0.55 Ma, with the Lau Biogeochemical Event starting within the hiatus (Arts et al., 2024). Our age for the onset of the Lau biogeochemical Event (424.06 ± 0.55 Ma) falls within this range, supporting alignment with previous results. The uncertainties assigned to durations in the record are around ~15%. Since only one proxy was used, a multi-proxy study might even further reduce the uncertainty in the underlying astrochronology of the Kosov quarry section.
The CWT in the time domain contains spectral peaks which can be linked to the 8.5-kyr half precession, 16-kyr precession, 28-kyr obliquity,50-kyr eccentricity, 100-kyr eccentricity, 200-kyr eccentricity, 405-kyr eccentricity and the 1.6-Myr eccentricity cycle (see Figure 7). The observation of a whole suite of spectral peaks corresponding to known astronomical cycles in the CWT scalogram validates the underlying age model. The periods of the obliquity and precession cycles are at the shorter end of the typically expected durations for the Silurian, yet remain within the current uncertainty range for the estimated durations of these cycles during the Silurian (Waltham, 2015; Farhat et al., 2022; Wu et al., 2024). The spectral peak of ~55 kyr was observed in the average CWT scalogram. This spectral peak can be either the 55 or 54-kyr eccentricity cycle or the 52.8-kyr obliquity cycle (Laskar et al., 2004; Laskar, 2020, p. 20) (see Figure 7). The spectral power of the 55-kyr cycle is highly variable, indicating that the 55-kyr cycle might be of a transient nature unrelated to astronomical forcing. A 200-kyr cycle can also be identified in the record. Given the strong imprint of the 100-kyr and 405-kyr eccentricity cycles, it is reasonable to interpret the 200-kyr cycle as the seldom observed 200-kyr eccentricity cycle (Hilgen et al. (2020)). The average spectral power of the induration record exhibits a peak at ~1.6 million years. This peak is most likely related to the 1.6-Myr eccentricity cycle. This peak is most likely related the 1.6-Myr eccentricity cycle Hinnov (2000); Laskar et al. (2004). Since the record is quite short (1800 kyr), its interpretation remains preliminary.
Since amplitude modulating cycles apply unidirectional to the amplitude of the signal, irrespective of whether the proxy exhibits a positive or negative response to insolation forcing it is possible to infer the phase relationship between astronomical cycles extracted from the proxy record and its relationship with the true phase of the astronomical cycle, as well as its relationship with changes in insolation (Hinnov, 2000; Laskar et al., 2004). The 405-kyr eccentricity cycle directly extracted from the record and derived from the Hilbert transform of the 100-kyr cycle shows an anti-phased relationship. The 100-kyr eccentricity cycle obtained directly from the induration record and from the Hilbert transform of the precession cycle also demonstrates a generally anti-phased relationship. This indicates that during isolation minima, induration values are at a maximum. This phase relationship is logical because during an isolation maximum, seasonality increases, leading to higher runoff, increased detrital products, and decreased carbonate productivity (lower induration) (Mutterlose and Ruffell, 1999; Martinez, 2018).
The established relationship between isolation and induration can also be used to understand the imprint of astronomical cycles in the rate of change ‰/kyr) of the δ¹³Ccarb record and the lag-1 sea-level curve. Since the 100-kyr eccentricity cycle extracted from the rate of change (‰/kyr) record is generally anti-phased with the 100-kyr eccentricity cycle extracted from the induration record, this indicates that the rate of change (‰/kyr) is the highest during eccentricity maxima. The 405-kyr eccentricity cycle extracted from the rate lag-1 record is generally anti-phased with the 405-kyr eccentricity cycle extracted from the induration record, indicating that sea-level is the highest during eccentricity maxima.
4.2 Astronomical cycles pacing the Ludfordian carbon cycle
The new age model enables the quantification of the rate of change in the δ¹³Ccarb record, yielding maximum values of 0.09‰ per kyr for the Lau Event. These rates were derived from a LOWESS regression of the δ¹³Ccarb curve to minimise the influence of short-term local noise. As such, the values represent robust regional and potentially global benchmarks for the Lau event. The high rate of change in the δ¹³Ccarb record (0.09‰ per kyr) during the Lau Event suggests that the Silurian carbon cycle was highly dynamic during the Lau Event. Since it is now possible to constrain the rate of change, it will be possible to model the Lau biogeochemical event in the domain. This will enable us to examine whether models, such as those of Frýda et al. (2020) till hold or whether alternative models that allow for non-steady behaviour need to be used. The 100-kyr eccentricity imprint in the δ¹³Ccarb rate-of-change record during the Lau Event indicates that the carbon cycle was at least partly paced by this cycle. The phase relationship shows that the highest rates of change occurred during eccentricity maxima. This suggests that eccentricity maxima resulted in enhanced organic carbon burial during the Lau Event. The observed phase relationship aligns with a scenario in which eccentricity maxima increased the amplitude of climate precession, resulting in enhanced seasonality (Rossignol-Strick, 1985; Van Os et al., 1994; Sachs and Repeta, 1999; Wang, 2009; Gambacorta et al., 2018). The resulting increase in seasonality intensified the hydrological cycle, increasing physical and chemical weathering, supplying fine-grained sediments and nutrients through enhanced (monsoonal) rainfall. Increased fluvial discharge promoted seasonal productivity blooms, facilitating the accumulation of organic-rich deposits. In contrast, during eccentricity minima, the reduced precession amplitude weakened monsoon intensity, lowered nutrient delivery, and re-established more oxic conditions.
4.3 Astronomical pacing the lag-1 sea-level curve
The lag-1 curve used as a sea-level proxy aligns with previous sea-level interpretation made for the Lau Biogeochemical Event (Manda and Kříž, 2006; Frýda et al., 2020) (see figure 10). The Lau Event is associated with a large regression and a period of low sea-level (see figure 10). The Lau Event is associated with a large regression and a period of low sea-level (see figure 10). The lag-1 curve indicates that a large regression occurred during the Lau which is in-line with previous interpretations which inferred that the Lau Event resulted in the Mid-Ludfordian Glaciation (Frýda et al., 2020). The 405-kyr eccentricity cycle extracted from the rate lag-1 record is generally anti-phased with that from the induration record, indicating that sea level is highest during eccentricity maxima (Figure 10). The relationship between astronomical forcing and sea-level is much akin to that observed when sea-level was under a glacio-eustatic regime (Lourens and Hilgen, 1997; Li et al., 2018). The transition from relatively stable low-amplitude variability into pronounced eccentricity-paced oscillations during the Lau Event further underscores the coupling between astronomical forcing, glacio-eustatic change, and the global carbon cycle during the Lau Event. The increased amplitude of sea-level variability indicates that astronomically paced glacio-eustatic fluctuations were amplified by the climatic changes associated with the Lau Event, culminating in the Mid-Ludfordian Glaciation.
5. Conclusions
The new astrochronology from the Kosov Quarry provides the first high-resolution temporal framework for the Lau biogeochemical Event. Cyclostratigraphic analysis demonstrates that the Lau biogeochemical Event started at 424.06 ± 0.55 Ma and lasted 0.87 ± 0.03 Myr, while the associated LKB Event was much shorter at 30 ± 10 kyr. Astronomical forcing is strongly imprinted across multiple proxies. The induration record is imprinted by astronomical cycles ranging from the 9-kyr half-precession to the 405-kyr eccentricity cycle. The δ¹³Ccarb rate-of-change record contains a strong imprint of the 100-kyr eccentricity cycle, whereas the lag-1 sea-level curve displays a moderate imprint by the 405-kyr eccentricity cycle. The δ¹³Ccarb rate-of-change curve further documents maximum values of 0.09 ‰/kyr during the onset of the Lau, highlighting the highly dynamic character of the Silurian carbon cycle. The imprint of the 100-kyr eccentricity cycle in the δ¹³Ccarb rate-of-change also indicates that this cycle acted as a pacemaker of carbon cycle variability and carbon burial during the Lau Event. The strong imprint of the 405-kyr eccentricity cycle in the lag-1 sea-level and the inferred phase relationship indicate that sea-level was subjected to an astronomically paced glacio-eustatic regime. Together, the results indicate that astronomical forcing amplified feedbacks between climate, organic carbon burial, and glacio-eustasy, shaping the magnitude and duration of the Lau Event. More broadly, these findings suggest that eccentricity-paced mechanisms played a pivotal role is shaping the climate, sea-level and carbon cycle dynamics during the Silurian.
6. References
Arts, M., Corradini, C., Pondrelli, M., Pas, D., and Da Silva, A.C., 2024, Age and orbital forcing in the upper Silurian Cellon section (Carnic Alps, Austria) uncovered using the WaverideR R package: Frontiers in Earth Science, v. 12, p. 1–24, doi:10.3389/feart.2024.1357751.
Farhat, M., Auclair-Desrotour, P., Boué, G., and Laskar, J., 2022, The resonant tidal evolution of the Earth-Moon distance: Astronomy & Astrophysics, v. 665, p. L1.
Fatka, O., and Mergl, M., 2009, The “microcontinent”Perunica: Status and story 15 years after conception: Geological Society, London, Special Publications, v. 325, p. 65–101, doi:10.1144/SP325.4.
Frýda, J., Lehnert, O., Frýdová, B., Farkaš, J., and Kubajko, M., 2020, Carbon and sulfur cycling during the mid-Ludfordian anomaly and the linkage with the late Silurian Lau/Kozlowskii Bioevent: Palaeogeography, Palaeoclimatology, Palaeoecology, v. 564, p. 110152, doi:10.1016/j.palaeo.2020.110152.
Frýda, J., Lehnert, O., Joachimski, M.M., Männik, P., Kubajko, M., Mergl, M., Farkaš, J., and Frýdová, B., 2021, The Mid-Ludfordian (late Silurian) Glaciation: A link with global changes in ocean chemistry and ecosystem overturns: Earth-Science Reviews, v. 220, p. 103652, doi:10.1016/j.earscirev.2021.103652.
Frýda, J., and Manda, Š., 2013, A long-lasting steady period of isotopically heavy carbon in the late Silurian ocean: Evolution of the δ13 record and its significance for an integrated δ13, graptolite and conodont stratigraphy: Bulletin of Geosciences, p. 463–482, doi:10.3140/bull.geosci.1436.
Gambacorta, G., Menichetti, E., Trincianti, E., and Torricelli, S., 2018, Orbital control on cyclical primary productivity and benthic anoxia: Astronomical tuning of the Telychian Stage (Early Silurian): Palaeogeography, Palaeoclimatology, Palaeoecology, v. 495, p. 152–162, doi:10.1016/j.palaeo.2018.01.003.
Hilgen, F., Zeeden, C., and Laskar, J., 2020, Paleoclimate records reveal elusive 200-kyr eccentricity cycle for the first time: Global and Planetary Change, v. 194, p. 103296, doi:10.1016/j.gloplacha.2020.103296.
Hinnov, L.A., 2000, New Perspectives on Orbitally Forced Stratigraphy: Annual Review of Earth and Planetary Sciences, v. 28, p. 419–475, doi:10.1177/0741713604268894.
HORNÝ, R., 1955, Base vrstev kopaninských eß 1 na jihozápadním okraji vulkanické facie (Kosov u Berouna): Věstník Ústředního ústavu geologického, v. 30, p. 81–86.
Kříž, J., 1992, Silurian field excursions: Prague Basin (Barrandian), Bohemia: National Museum of Wales, v. 13.
Kříž, J., 1991, The Silurian of the Prague Basin (Bohemia)–tectonic, eustatic and volcanic controls on facies and faunal development: Special Papers in Palaeontology, v. 44, p. 179–203.
Laskar, J., Robutel, P., Joutel, F., Gastineau, M., Correia, a.C.M., and Levrard, B., 2004, A long-term numerical solution for the insolation quantities of the Earth: Astronomy and Astrophysics, v. 428, p. 261–285, doi:10.1051/0004-6361:20041335.
Lehnert, O., Frýda, J., Buggisch, W., Munnecke, A., Nützel, A., Křiž, J., and Manda, S., 2007, δ13C records across the late Silurian Lau event: New data from middle palaeo-latitudes of northern peri-Gondwana (Prague Basin, Czech Republic): Palaeogeography, Palaeoclimatology, Palaeoecology, v. 245, p. 227–244, doi:10.1016/j.palaeo.2006.02.022.
Li, M., Hinnov, L.A., Huang, C., and Ogg, J.G., 2018, Sedimentary noise and sea levels linked to land-ocean water exchange and obliquity forcing: Nature Communications, v. 9, doi:10.1038/s41467-018-03454-y.
Lourens, L.J., and Hilgen, F.J., 1997, Long-periodic variations in the earth’s obliquity and their relation to third-order eustatic cycles and late Neogene glaciations: Quaternary International, v. 40, p. 43–52, doi:10.1016/S1040-6182(96)00060-2.
Manda, Š., and Kříž, J., 2006, Environmental and biotic changes in subtropical isolated carbonate platforms during the Late Silurian Kozlowskii Event, Prague Basin: GFF, v. 128, p. 161–168.
Martinez, M., 2018, Mechanisms of Preservation of the Eccentricity and Longer-term Milankovitch Cycles in Detrital Supply and Carbonate Production in Hemipelagic Marl-Limestone Alternations: Elsevier Ltd, v. 3, 189–218 p., doi:10.1016/bs.sats.2018.08.002.
Meyers, S.R., 2019, Cyclostratigraphy and the problem of astrochronologic testing: Earth-Science Reviews, v. 190, p. 190–223, doi:10.1016/j.earscirev.2018.11.015.
Meyers, S.R., 2015, The evaluation of eccentricity-related amplitude modulation and bundling in paleoclimate data: An inverse approach for astrochronologic testing and time scale optimization: Paleoceanography, v. 30, p. 1625–1640, doi:10.1002/2015PA002850.
Mutterlose, J., and Ruffell, A., 1999, Milankovitch-scale palaeoclimate changes in pale-dark bedding rhythms from the Early Cretaceous (Hauterivian and Barremian) of eastern England and northern Germany: Palaeogeography, Palaeoclimatology, Palaeoecology, v. 154, p. 133–160, doi:10.1016/S0031-0182(99)00107-8.
Sachs, J.P., and Repeta, D.J., 1999, Oligotrophy and nitrogen fixation during eastern Mediterranean sapropel events: Science, v. 286, p. 2485–2488.
Scotese, C.R., 2021, An atlas of phanerozoic paleogeographic maps: The seas come in and the seas go out: Annual Review of Earth and Planetary Sciences, v. 49, p. 679–728, doi:10.1146/annurev-earth-081320-064052.
Štorch, P., 1995, Biotic crises and post-crisis recoveries recorded by Silurian planktonic graptolite faunas of the Barrandian area (Czech Republic): Geolines, v. 3, p. 59–70.
Tasáryová, Z., Schnabl, P., Čížková, K., Pruner, P., Janoušek, V., Rapprich, V., Štorch, P., Manda, Š., Frýda, J., and Trubač, J., 2014, Gorstian palaeoposition and geotectonic setting of Suchomasty Volcanic Centre (Silurian, Prague Basin, Teplá-Barrandian Unit, Bohemian Massif): Gff, v. 136, p. 262–265, doi:10.1080/11035897.2013.879735.
Van Os, B., Lourens, L., Hilgen, F., De Lange, G., and Beaufort, L., 1994, The formation of Pliocene sapropels and carbonate cycles in the Mediterranean: Diagenesis, dilution, and productivity: Paleoceanography, v. 9, p. 601–617.
Waltham, D., 2015, Milankovitch Period Uncertainties and Their Impact On Cyclostratigraphy: Journal of Sedimentary Research, v. 85, p. 990–998, doi:10.2110/jsr.2015.66.
Wang, P., 2009, Global monsoon in a geological perspective: Chinese Science Bulletin, v. 54, p. 1113–1136.
Wu, Y., Malinverno, A., Meyers, S.R., and Hinnov, L.A., 2024, A 650-Myr history of Earth’s axial precession frequency and the evolution of the Earth-Moon system derived from cyclostratigraphy: Science Advances, v. 10, p. eado2412.
